## setup ##
  
  # clean environment
  rm(list=ls())
  
  # load packages
  library(xlsx)
  library(sandwich)
  library(lmtest)
  library(tidyverse)
  library(xtable)
  library(readxl)
  library(stargazer)
  
  # load main data
dat <- read.csv2("data_experiment_w_late_replies_csv.csv",
                 encoding = 'UTF-8')
  

## assign correct types (just to be safe)
# outcome variable to numeric
dat$response <- as.numeric(as.character(dat$response))

# IVs to numeric
dat$immback <- as.numeric(as.character(dat$immback))
dat$partisan <- as.numeric(as.character(dat$partisan))


## Figure A2: Main Regression with different day samples

ests <- tibble()
limits <- c(14, 21, 28, 35, 42, 363)

for (i in 1:length(limits)) {
  print(limits[i])
  dat_loop <- dat %>% 
    mutate(response = ifelse(day_diff_notimputed > limits[i] | is.na(day_diff_notimputed) , 0, response )) # if day_diff_notimputed is bigger than i, set 0

# Model without covariates:
   f1_overall <- lm(response ~ immback + partisan , dat_loop)
  
  # robust standard errors
  f1_overall_rob <- broom::tidy(coeftest(f1_overall, vcov = vcovHC(f1_overall, "HC1"))) 
  
  f1_overall_rob <- f1_overall_rob %>% 
    filter(term %in% c("immback","partisan")) %>% 
    mutate(cutoff = paste(limits[i], "days"),
           conf.low95 = estimate -  qnorm(0.975) * std.error,
           conf.high95 = estimate + qnorm(0.975)* std.error)
  
  ests <- rbind(ests, f1_overall_rob )
}

ests <- ests %>% drop_na(term) 
ests <- ests %>%    mutate(cutoff = ifelse(cutoff == "363 days", ">50 days", cutoff)) %>% 
  mutate(cutoff = as.factor(cutoff))

ests$cutoff <- ordered(ests$cutoff, levels = c("14 days",
                                              "21 days",
                                              "28 days",
                                              "35 days",
                                              "42 days",
                                              ">50 days"))

pd <- position_dodge(1)

ests <- ests %>% 
  mutate(term = dplyr::recode(term,
                              `partisan` = 'Partisan',
                              `immback` = 'Immigrant Background'))  %>% 
  rename(Coefficient = estimate,
         Period = cutoff)

p1 <- ggplot(ests, aes(Period, Coefficient)) +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  geom_errorbar(aes(ymin = conf.low95, ymax = conf.high95),
                # position = pd,
                width=.2, size=.4) +
  geom_point() +
  facet_wrap(~term, scales = 'free_y', ncol = 2) +
  theme_bw() +
  xlab("Observation Period") + ylab("Treatment Effect") +
  theme(panel.grid = element_blank(), legend.position = "bottom", text = element_text(size=14)) 

p1

# Uncomment to save plot as PDF:
# ggsave('fig_A2.pdf',
#        width=9, height=5)

